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Abstract 

We propose a hybrid dynamical system approach to model the evolu- 
tion of a pathogen that experiences different selective pressures according to 
a stochastic process. In every environment, the evolution of the pathogen 
is described by a version of the Fisher-Haldane- Wright equation while the 
switching between environments follows a Markov jump process. We investi- 
gate how the qualitative behavior of a simple single-host deterministic system 
changes when the stochastic switching process is added. In particular, we 
study the stability in probability of monomorphic equilibria. We prove that 
in a "constantly" fluctuating environment, the genotype with the highest 
mean fitness is asymptotically stable in probability while all others are un- 
stable in probability. However, if the probability of host switching depends 
on the genotype composition of the population, polymorphism can be stably 
maintained. 

Remark. This is a corrected version of the paper that appeared in 
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1. Introduction 

Living organisms face changing environmental conditions. Parasites are 
a case in point: after each transmission event they find themselves in a new 
host that may be quite different from the previous one. For example, the 
immune system of the new host may respond differently to the parasite, and 
the new host may have a different genotype or even belong to a different 
species. Thus, a parasite genotype that is characterized by a high fitness in 
one host may have a low fitness in a different host. The question therefore 
arises how parasites evolve under the fluctuating selective pressures imposed 
on them through transmission events to different hosts. 

Most of the studies so far have focused on models for host-pathogen in- 
teractions in a deterministic context [H1II2]- In some applications however it 
is natural to assume that environment (and hence fitness landscape) switch- 
ing is not deterministic. For example, a pathogen could switch to a different 
host. Evolution of the pathogen then takes place in the new host (or environ- 
ment), where the pathogen genotypes face different selective pressures, hence 
the dynamics of the pathogen genotypes are different. We remark that the 
evolving organism need not be a pathogen, nor is the environment necessarily 
a "host". 

Evolution of organisms in deterministically and randomly varying envi- 
ronments has been studied by many authors, see [2] for an early review. 
Karlin and collaborators [5J E] introduced both deterministic and stochastic 
models for the evolution of haploid and diploid organisms under changing 
selection intensities for fixed and varying population sizes. In case of a deter- 
ministic two-allele model they showed that the genotype with higher selection 
intensity goes to fixation and the time to fixation varies according to the se- 
lection intensities. Furthermore, they investigated a stochastic model where 
generational selection intensities are identically distributed independent ran- 
dom variables. They focused on the question how the probabilities of fixation 
and the times to fixation change in the stochastic model. In |7j, Kirzhner et 
al. considered a 4-dimensional system of difference equations for the haplo- 
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type frequencies of a two-locus model. Typical two-locus models show either 
fixation in one or both loci or stable polymorphic cycles, with period equal- 
ing the period of the environmental changes, i.e. the periodic fitness values. 
They however showed the existence of so called supercycles that have 1100 
times the period of the periodically changing environment. The questions 
of structural stability, i.e. sensitivity in terms of the fitness parameters and 
the size of the basin of attraction of these cycles were investigated. Simi- 
larly, Nagylaki [10J investigated the existence of genetic polymorphisms for 
two-allele models with periodically varying fitness values. He showed that in 
a continuous differential equation model genetic polymorphism will persist 
with periods equaling the periods of the varying fitness values, however in a 
discrete model fixation is also possible. 

Hybrid switching differential equations and more generally hybrid switch- 
ing diffusions have found many applications in wireless communications, 
queuing networks, ecology (15] and financial mathematics, to name but a few; 
see [2] and the references therein. The word "hybrid" refers to the coexis- 
tence of continuous dynamics and discrete events, see also the related concept 
of piecewise deterministic processes pp. In this paper we study a simplified 
version of the continuous time Fisher-Haldane- Wright equation (also known 
as standard replicator equation, [3]) subject to fitnesses driven by a Markov 
jump process. It is well known that in the deterministic Fisher-Haldane- 
Wright model the pathogen genotype that has the highest fitness value will 
go to fixation. The coupling of different Fisher-Haldane- Wright equations by 
a Markov process however requires a new definition of the concept of "highest 
fitness" . Hence we study the possible changes to the stability behaviors of the 
monomorphic equilibrium states depending on the stationary distribution of 
the switching Markov process. First we establish analytical results for the 
stability/instability of equilibria in the hybrid model. We show that in the 
case of a state-independent switching process, the monomorphic equilibrium 
of the genotype with highest mean fitness is asymptotically stable in proba- 
bility while the monomorphic equilibria of all other genotypes are unstable in 
probability. As the stationary distribution of the switching Markov process 
varies, so does the mean fitness of each genotype. This results in exchanges 
of stability without the merger of the equilibria during the transition process. 
We may call this a "stochastic transcritical bifurcation" . Finally, we present 
some numerical simulations to illustrate our result and also an example of a 
state-dependent switching process. 
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2. The switching differential equation 



We consider a model for m genotypes of a pathogen evolving in n pos- 
sible environments. Let w% > denote the fitness value of genotype % in 
environment k. We assume for simplicity that for any fixed environment k, 
the fitness values wf are all different. We write Wj for the vector of all fitness 
values of genotype i. Let Pi denote the frequency of pathogen i, so that the 
dynamics in each environment takes place in the (m— l)-dimensional simplex 



i=i 



We write P(t) = (Pi(t), P2 (t), ■ ■ ■ , P m (t)) for the state of the system at time 
t. The frequency dynamics of pathogen genotype i in environment k is given 
by 

dPi(t) 



dt 



Piit) ( w{ - Y.wfcit) ) =: F t (P(t),k). (2.1) 



This equation is the Fisher-Haldane- Wright equation for frequencies of geno- 
types of asexually proliferating organisms. The rate of growth or decay of a 
genotype is determined by the difference of its fitness and the average fitness 
of the population. Observe that the simplex T m_1 and any of its subsimplices 



are invariant under the dynamics given by equation (2.1). It can be shown 



by straightforward computation that the average fitness in environment k 

m 

fe (P) = 5>* P * ( 2 - 2 ) 



i=l 



satisfies 
d 



r (p(t)) = J2 w i p i w * -J2 w J p j = E p *K fc ) 2 - E p < 
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with equality if and only if P is an equilibrium. It follows from the global 



existence of solutions and LaSalle's theorem that every trajectory of (2.1) 
approaches one of the finitely many equilibria situated at the vertices of the 
simplex, see [8]. 

The environment switches according to a continuous time stochastic pro- 
cess a(t) that takes values in the set M. — {1, 2, . . . , n}. The switching pro- 
cess a is a Markov process with (possibly state-dependent) generator matrix 
Q(P) whose entries qki(P) are defined by 

P{a(t+At) = 1 1 a(t) = k, (Pis), a(s)), s<t} = q kl (P(t))At+o(At). (2.3) 

The elements q k i of the generator matrix Q satisfy q^i > for all k ^ I and 
q^i = for every k G M. (such a matrix is said to have the q-property, 

leM 

see [13J). The complete hybrid switching ordinary differential equation can 
be cast in the form 

dP 

*=w>. «<«)). (24) 

P(0) = P eT m -\ a(0) = aeM, a. s., 
where a(t) = k determines the environment k at time t and F = (iq, . . . , F m ) 



are defined by (2.1). The right hand side of the differential equation in 
(2.4) is globally Lipschitz continuous on the compact set T m_1 x Ai. This 
implies global existence and uniqueness of solutions in the sense of stochastic 
processes, see [HI Theorem 2.1]. 



For the equation (2.1) restricted to a fixed environment k, the vertices 
ei of the simplex T m_1 (i.e. the unit vectors of lR m ) are all the equilibrium 
solutions. It is easy to show that all but one of these equilibria are unstable 
and that the stable equilibrium in environment k is the one for which the 
fitness value wf is the largest. In the following section we investigate how 
this result generalizes to the case that stochastic switching is introduced. For 
this of course, we need to first generalize the concept of stability to switching 
ordinary differential equations. 



3. Stability and instability in probability 

In this section we establish results concerning the stability and instabil- 
ity of monomorphic steady states of the hybrid model. We first recall the 
following definition [HI Definition 8.1]. 
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Definition 3.1. Let (X x ' a (t)) t > be the solution of a hybrid switching ordi- 
nary differential equation 



X(t) = F(X(t),a(t)), 
P{a(t + At) = l\ a(t) = k, (X(s), a(s)), s < t} = q M (X(t))At + o{At), 
X(0) = x, a(0) = a a. s., 



and let (without loss of generality) x = be an equilibrium solution, i.e. a 
solution of the equation F(0,a) = for every a £ Ai. We say that is 
stable in probability if 



for every a £ Ai and every r > 0. We say that is asymptotically stable 
in probability if it is stable in probability and 



for every a £ Ai. Finally, is unstable in probability if it is not stable 
in probability. 

For n-tuples of functions g( ■ , k) £ C 1 (IR m ) one defines a linear operator 
£, the stochastic Lie derivative (see [HI Equation (8.3), p. 219]) 



where V denotes the gradient with respect to the x- variable for fixed k £ Ai. 
This is a natural generalization of the derivative of a scalar function along a 
vector field well known in the theory of ordinary differential equations. The 
following is Proposition 8.6, [HJ p. 223]. 

Theorem 3.2. Let D C IR m be a neighborhood of and assume that there 
exists a function V : D x Ai — > [0, oo) with the following properties 

• V( • , k) is continuous and vanishes only at 0, 

• V( • , k) is continuously differentiable in D\ {0}, and 





n 



Cg(x, k) = F(x, k) ■ Vg(x, k) + ^ Qki(x)g(x, k), 



i=i 
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• there exists a function k : (0,r) — > (0, oo) such that for all k G Ai and 

\x\ > g, 

CV(x,k) < -k(q) < 0. 
Then the equilibrium x = is asymptotically stable in probability. 

A function that satisfies the conditions of the theorem is called a Lyapunov 
function (for asymptotic stability). 

Throughout the remainder of this section we consider the case of a state- 
independent generator matrix Q with a universal stationary distribution 7r = 
(qi, . . . , q n ). This is the solution of the equations 

7i ■ 1 = 1, and ttQ = 0. 

If qki > for all k I then the matrix Q is irreducible and tt > is unique 

pa p. 2i]. 



Theorem 3.3. Let P\ be the genotype with the highest mean fitness, that is 
7r ■ wi > 7r ■ Wj for all i = 2, . . . , m (3.5) 
Then the equilibrium e\ is asymptotically stable in probability. 

Remark 3.4 For almost every stationary distribution ix £ T™" 1 , exactly 



one genotype satisfies a condition similar to (3.5). 

Proof. For i — 2, . . . , m we set a*j = wf — w\ for the difference of fit- 
ness values with respect to genotype 1 and a^i = (a} l7 . . . , a"i). Using the 

m 

constraint Pj (t) = 1, we eliminate Pi and obtain the reduced systems 



dPi(t) _ k 



for i = 2, . . . , m and k = 1, . . . , n. Notice that for fixed environment k the 
linear part of this system has a diagonal structure. We define 

Pi := -7r • a^i > 0, 
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with the last inequality holding true since genotype 1 has the higher mean 
fitness compared to every other genotype. For i = 2, ... ,m we solve the 
systems of equations 

Qci = a^i 

for the vector Cj = (cj, . . . , cf) where 1 is the column vector with n entries 
1. The right hand sides of these equation are orthogonal to the kernel of Q 
which is spanned by 1, hence there exist solutions. For i = 2, . . . , m and 
k — 1, . . . , n, we define 

V l {P l1 k) = (l- 1 c k i )Pl P t >0, 

with < 7 < 1 yet to be selected, in such a way that all coefficients are 
positive. We have 



w(p t , k) = 7 (i - i$)pr\<iPi + o(i)) + Yl mi - vbp* 



= 7*7 f (1 - 14)4,1 ~ E 9*4 + ( 3 - 6 ) 

= 7^ 7 ((1 - ic^al - (al + A) + o(l)) 
= 7 ^ 7 (-7C- aji + 7T ■ a,,! + o(l)) , 

where we have made use of the fact that the row sums of Q are zero. In order 
to make all the factors in parentheses negative, we have to choose < 7 < 1 
such that the inequality 

7T ■ a M < 7cfaf (3.7) 



holds. By assumption (3.5), the left hand side of inequality (3.7) is negative. 
Therefore, for those indices i and k for which cfa^ > 0, no condition arises 
for 7. If on the other hand cfaf 1 < 0, then we can select 

h a i,i 

Although the cf are not explicitly known, this is a minimum of finitely many 
positive numbers. The Lyapunov function is the sum of functions of a single 
variable 



V{P 2 ,...,P m ,k) = Y i V i {P i ,k) 



i=2 



and the condition of Proposition 8.6 in [13] follows from the linearity of the 
operator C and the choice of 7. □ 
Instability in probability of an equilibrium can be proved similarly The 
following is Proposition 8.7, [HI p. 223]. Notice however that the Lyapunov 
function does not vanish but has a pole at the unstable equilibrium. 

Theorem 3.5. Let D C IR m be a neighborhood of and assume that there 
exists a function V : D x M. — > [0, 00) with the following properties 

• V{ ■ , k) is continuously differentiable in D\ {0}, and 

• there exists a function k : (0,r) — > (0, 00) such that for all k E M and 

\x\ > g, 

CV(x,k) < -k(q) < 0, 

• for all k e M, 

lim V(x, k) = 00. 

|x|->0 

Then the equilibrium x = is unstable in probability. 



Theorem 3.6. Under the assumption (3.5) of Theorem 3.3, the equilibrium 
ej, i — 2, . . . , m is unstable in probability. 



Proof. The proof is very similar to that of Theorem |3.3[ so we only give 
a sketch here. This time it is Pi that is being eliminated from the system 
containing Pi and Pj. This results in the reduced systems 

for / 7^ i and a\ i = wf — wf. For i = 2, . . . , m let Cj = (cj, . . . , c") be the 
solution of 

Qci = ai ; j — PjI. 

We set 

V{P 1 , . . . , P^, P i+l , . . . , P m , k) = V(P 1 , k) = (1 - ^Dpi P 1 > 0, 
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where > 7 > — 1 has yet to be selected, small enough that all coefficients 
are positive. With a calculation similar to (3.6) we obtain 



CV{Pi, k) = 7 (1 - 7 cf)Pr 1 «,Pi + o(l)) + QhiO- - Tt)Pi 



3=1 



= l P i ((! _ l c i) a i,i - ( a i,i _ + 
= 7 P7 (-7C- ^ + vr • a M + o(l)) . 

In order to make all the factors in parentheses positive (so that the entire 
expression becomes negative), we need to have 

n J 71 ' a M k k n 1 

> 7 > max <^ : <0>- 



i—2,. ..m 
fc — 1 n 



era 



The expressions whose maximum is taken are all negative since it ■ &i , > by 



assumption (3.5). The condition of Proposition 8.7 in [13] is thereby verified. 
□ 

Remark 3.7 The notion of "highest mean fitness" requires that the gener- 
ator matrix Q is independent of the state and so has a universal stationary 
distribution n. If Q depends continuously on P, it is still possible to formu- 
late the corresponding "local" stability results for the equilibria by taking 
7r to be a stationary distribution of Qfa). 

4. Numerical simulations and examples 

The following is an interesting example of how stability can arise through 
stochastic coupling. Let the fitness values of three genotypes in two environ- 
ments be given by 

1 ! 1 7 ! 11 

W\ = 1, Wn = , Wn = , 

1 ' 2 10' 3 10' 

2 1 2 H 2 7 
Wi = 1, U), = , II) = . 

1 ' 2 10' 3 10 
Although genotype 1 does not have the highest fitness in any environment, 
it has the highest mean fitness for stationary distributions (q, 1 — q) with 
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I < q < |, see Figure [TJ If the generator matrix Q of the Markov process does 
not depend explicitly on the state P(t), we can determine the switching times 
a priori according to ti+i = U + r, where r is an exponentially distributed 
random variable with mean 1 (for example). 





Figure 1: (Left panel) The mean fitness of the three equilibria as a function of the pa- 
rameter q of the stationary distribution. (Right panel) Trajectories converging to e\ when 
q = i Parts of the trajectory in red indicate that environment 1 is active (when e% is at- 
tracting) while parts of the trajectory in blue indicate that environment 2 is active (when 
e2 is attracting). 



To finish this section, we present two example with a state-dependent 
generator matrix Q(P). Let n = m = 2, 



and define two switching matrix functions 



10' 





Pi 


Pi 


-Pi 


-Pi 


Pi 


Pi 


-Pi 



QiiPuPi) 



Q2(Pl,P 2 ) 



This choice of the generator matrix means that the jump process favors jumps 
into the environment that is beneficial (in case Qi), respectively disadvan- 
tageous (in case Q 2 ) for the genotype that currently dominates. In contrast 
to the previous simulations with state-independent generator matrix, it is 
now necessary to update the transition matrix of the Markov chain, namely 
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exp(Q(P)At) during each time step of length At. Following [TJ1 Chapter 
5.3], we use the approximation / + Q(P)At. The stationary distribution of 
Qi(ei) is, incidentally, for i = 1,2. It follows from Theorem 3.3 and Re- 



mark [3^7] that both equilibria are locally asymptotically stable in probability. 
Conversely, the stationary distribution of Q2(ei) is €2 and vice versa. Under 
this regime, both equilibria are locally unstable in probability. The results 
in Figure [2] show that stochastic bistability may arise (for the choice Q\(P), 
left panel) or that solutions do not converge to a monomorphic steady state 
(for the choice Q2(P), right panel). 




Figure 2: A state-dependent generator matrix Qi(P) leads to bistability (left panel) 
whereas matrix Q2{P) results in failure to converge to an equilibrium (right panel). Parts 
of the trajectory in red indicate that environment 1 is active (when e\ is attracting) while 
parts of the trajectory in blue indicate that environment 2 is active (when is attracting). 



In terms of biological interpretation, one can conceive of competing pathogen 
genotypes that cause different behaviors in the affected host. For example, if 
the dominating pathogen genotype has only mild effects on their host's well- 
being, infected individuals may retain their usual mobility and thereby make 
a transition into a new environment more likely. On the other hand, if the 
dominating pathogen genotype causes severe morbidity, the host may exhibit 
restricted mobility so that a transition into a new environment becomes less 
likely. 
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5. Conclusions 



In this work we consider the dynamics of a simple host-pathogen system, 
where pathogen genotype frequencies evolve according to a simple determin- 
istic model. The selective pressures switch according to a Markov process. 
We use the framework of switching differential equations to compare the evo- 
lution of the pathogen in a single deterministic versus a hybrid system. In 
the switching system interesting new stability patterns emerge, depending 
on the stationary distribution of the underlying Markov process. We assume 
a fixed number of environments (and corresponding fitnesses), in contrast to 
previous works. For example, Karlin and collaborators (5J E] assumed that 
the fitnesses during each generation are independent identically distributed 
random variables. Gillespie on the other hand in jl] proposed a stochastic 
differential equation where the fitness is a process with continuous sample 
paths. 

In the case of a state-independent generator matrix of the Markov pro- 
cess Q, we have a partition of the simplex T n_1 of all possible stationary 
distributions tc into regions where one genotype has a greater mean fitness 
than all others, except for a set of measure zero where two genotypes have 
equal mean fitness (where bifurcations occur). This complete classification 
relies on the diagonal structure of the Jacobian of the reduced system (??) 
at the equilibrium 0. Due to this decoupling, it is possible to use a sum 
of Lyapunov functions that all depend on one variable only. In this way, 
we obtain a condition for asymptotic stability using convex combinations of 
corresponding elements of the spectra of the Jacobians in the different envi- 
ronments. We expect such a result to hold in the greater context of switching 
ordinary differential equations and diffusion processes with regime switching. 

Our work can be refined and extended in various ways. Firstly, we use a 
very simple deterministic competition model (2.1), where the pathogen geno- 
types are ordered according to their fitness values and the only equilibria are 
the vertices of the simplex T™^ 1 . A straightforward extension would be 
to consider the continuous time Fisher-Haldane- Wright equation for diploid 
organisms for which there exist equilibria in the interior of T m_1 . Other 
competition models may lead to deterministic bistability or to periodic or- 
bits (for example the rock-paper-scissors game [3]). Secondly, although the 
host switching process is stochastic, we model within-host evolution in a de- 
terministic way. A more realistic approach would incorporate random genetic 
drift into the model. This may be particularly important during transmis- 
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sion events, which often involve population bottlenecks due to small inoculum 
sizes. Finally, our model only considers a single chain of transmission events 
and neglects between-host selection as well as superinfections. It may be 
possible to also consider multiple (branching and coalescing) transmission 
chains and thus fully couple within-host and epidemiological dynamics [Sj. 
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